home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / viewers / polyview / polyvw31.lha / Polyview3.1 / new / quantizer.c < prev    next >
C/C++ Source or Header  |  1993-06-23  |  15KB  |  479 lines
  1. /*****************************************************************************
  2.  * NCSA Polyview 3.0                                                         *
  3.  *                                                                           *
  4.  * Version 3 changes and additions by Marc Andreessen.                       *
  5.  * Version 2 by Brian Calvert.                                               *
  6.  *                                                                           *
  7.  * Software Development Group                                                *
  8.  * National Center for Supercomputing Applications                           *
  9.  * University of Illinois at Urbana-Champaign                                *
  10.  *                                                                           *
  11.  * This is BETA release software.  As such it may contain software bugs and  *
  12.  * exhibit inconsistencies.                                                  *
  13.  *                                                                           *
  14.  * Please send bug reports to polyview@ncsa.uiuc.edu.                        *
  15.  *                                                                           *
  16.  * Copyright (c) 1992 The Board of Trustees of the University of Illinois.   *
  17.  *                                                                           *
  18.  * Permission to use, copy, and modify this software and its                 *
  19.  * documentation for educational, research, and non-profit purposes is       *
  20.  * hereby granted, provided that the above copyright notice, the original    *
  21.  * authors names, and this permission notice appear in all such copies.      *
  22.  * Any distribution of this software requires the explicit and written       *
  23.  * authorization of the authors.                                             *
  24.  *                                                                           *
  25.  * The University of Illinois makes no representations about the             *
  26.  * suitability of this software for any purpose.  It is provided "as is"     *
  27.  * without warranty of any kind.                                             *
  28.  *****************************************************************************/
  29.  
  30. /* $Header: /usr3/people/gbourhis/pv3/new/RCS/quantizer.c,v 1.2 92/09/19 06:30:59 marca Exp $ */
  31.  
  32. #ifdef RCSLOG
  33. $Log:    quantizer.c,v $
  34.  * Revision 1.2  92/09/19  06:30:59  marca
  35.  * Added bzeros for repeated use.
  36.  * 
  37.  * Revision 1.1  1992/09/18  10:55:26  marca
  38.  * Initial revision
  39.  *
  40. #endif
  41.  
  42. /* Our sincere thanks to the author of this file (see notes below) for
  43.    making his code available for redistribution. */
  44.  
  45. /**********************************************************************
  46.         C Implementation of Wu's Color Quantizer (v. 2)
  47.         (see Graphics Gems vol. II, pp. 126-133)
  48.  
  49. Author:    Xiaolin Wu
  50.     Dept. of Computer Science
  51.     Univ. of Western Ontario
  52.     London, Ontario N6A 5B7
  53.     wu@csd.uwo.ca
  54.  
  55. Algorithm: Greedy orthogonal bipartition of RGB space for variance
  56.        minimization aided by inclusion-exclusion tricks.
  57.        For speed no nearest neighbor search is done. Slightly
  58.        better performance can be expected by more sophisticated
  59.        but more expensive versions.
  60.  
  61. The author thanks Tom Lane at Tom_Lane@G.GP.CS.CMU.EDU for much of
  62. additional documentation and a cure to a previous bug.
  63.  
  64. Free to distribute, comments and suggestions are appreciated.
  65. **********************************************************************/    
  66.  
  67. #include <stdio.h>
  68.  
  69. #define MAXCOLOR    256
  70. #define    RED    2
  71. #define    GREEN    1
  72. #define BLUE    0
  73.  
  74. struct box {
  75.     int r0;             /* min value, exclusive */
  76.     int r1;             /* max value, inclusive */
  77.     int g0;  
  78.     int g1;  
  79.     int b0;  
  80.     int b1;
  81.     int vol;
  82. };
  83.  
  84. /* Histogram is in elements 1..HISTSIZE along each axis,
  85.  * element 0 is for base or marginal value
  86.  * NB: these must start out 0!
  87.  */
  88.  
  89. float        m2[33][33][33];
  90. long int    wt[33][33][33], mr[33][33][33],    mg[33][33][33],    mb[33][33][33];
  91. unsigned char   *Ir, *Ig, *Ib;
  92. int            quant_size; /*image size*/
  93. int        K;    /*color look-up table size*/
  94. unsigned short int *Qadd;
  95. unsigned char    lut_r[MAXCOLOR], lut_g[MAXCOLOR], lut_b[MAXCOLOR];
  96.  
  97. void
  98. Hist3d(vwt, vmr, vmg, vmb, m2) 
  99. /* build 3-D color histogram of counts, r/g/b, c^2 */
  100. long int *vwt, *vmr, *vmg, *vmb;
  101. float    *m2;
  102. {
  103. register int ind, r, g, b;
  104. int         inr, ing, inb, table[256];
  105. register long int i;
  106.         
  107.     for(i=0; i<256; ++i) table[i]=i*i;
  108.     Qadd = (unsigned short int *)malloc(sizeof(short int)*quant_size);
  109.     for(i=0; i<quant_size; ++i){
  110.         r = Ir[i]; g = Ig[i]; b = Ib[i];
  111.         inr=(r>>3)+1; 
  112.         ing=(g>>3)+1; 
  113.         inb=(b>>3)+1; 
  114.         Qadd[i]=ind=(inr<<10)+(inr<<6)+inr+(ing<<5)+ing+inb;
  115.         /*[inr][ing][inb]*/
  116.         ++vwt[ind];
  117.         vmr[ind] += r;
  118.         vmg[ind] += g;
  119.         vmb[ind] += b;
  120.              m2[ind] += (float)(table[r]+table[g]+table[b]);
  121.     }
  122. }
  123.  
  124. /* At conclusion of the histogram step, we can interpret
  125.  *   wt[r][g][b] = sum over voxel of P(c)
  126.  *   mr[r][g][b] = sum over voxel of r*P(c)  ,  similarly for mg, mb
  127.  *   m2[r][g][b] = sum over voxel of c^2*P(c)
  128.  * Actually each of these should be divided by 'size' to give the usual
  129.  * interpretation of P() as ranging from 0 to 1, but we needn't do that here.
  130.  */
  131.  
  132. /* We now convert histogram into moments so that we can rapidly calculate
  133.  * the sums of the above quantities over any desired box.
  134.  */
  135.  
  136.  
  137. void
  138. M3d(vwt, vmr, vmg, vmb, m2) /* compute cumulative moments. */
  139. long int *vwt, *vmr, *vmg, *vmb;
  140. float    *m2;
  141. {
  142. register unsigned short int ind1, ind2;
  143. register unsigned char i, r, g, b;
  144. long int line, line_r, line_g, line_b,
  145.      area[33], area_r[33], area_g[33], area_b[33];
  146. float    line2, area2[33];
  147.  
  148.     for(r=1; r<=32; ++r){
  149.     for(i=0; i<=32; ++i) 
  150.         area2[i]=area[i]=area_r[i]=area_g[i]=area_b[i]=0;
  151.     for(g=1; g<=32; ++g){
  152.         line2 = line = line_r = line_g = line_b = 0;
  153.         for(b=1; b<=32; ++b){
  154.         ind1 = (r<<10) + (r<<6) + r + (g<<5) + g + b; /* [r][g][b] */
  155.         line += vwt[ind1];
  156.         line_r += vmr[ind1]; 
  157.         line_g += vmg[ind1]; 
  158.         line_b += vmb[ind1];
  159.         line2 += m2[ind1];
  160.         area[b] += line;
  161.         area_r[b] += line_r;
  162.         area_g[b] += line_g;
  163.         area_b[b] += line_b;
  164.         area2[b] += line2;
  165.         ind2 = ind1 - 1089; /* [r-1][g][b] */
  166.         vwt[ind1] = vwt[ind2] + area[b];
  167.         vmr[ind1] = vmr[ind2] + area_r[b];
  168.         vmg[ind1] = vmg[ind2] + area_g[b];
  169.         vmb[ind1] = vmb[ind2] + area_b[b];
  170.         m2[ind1] = m2[ind2] + area2[b];
  171.         }
  172.     }
  173.     }
  174. }
  175.  
  176.  
  177. long int Vol(cube, mmt) 
  178. /* Compute sum over a box of any given statistic */
  179. struct box *cube;
  180. long int mmt[33][33][33];
  181. {
  182.     return( mmt[cube->r1][cube->g1][cube->b1] 
  183.        -mmt[cube->r1][cube->g1][cube->b0]
  184.        -mmt[cube->r1][cube->g0][cube->b1]
  185.        +mmt[cube->r1][cube->g0][cube->b0]
  186.        -mmt[cube->r0][cube->g1][cube->b1]
  187.        +mmt[cube->r0][cube->g1][cube->b0]
  188.        +mmt[cube->r0][cube->g0][cube->b1]
  189.        -mmt[cube->r0][cube->g0][cube->b0] );
  190. }
  191.  
  192. /* The next two routines allow a slightly more efficient calculation
  193.  * of Vol() for a proposed subbox of a given box.  The sum of Top()
  194.  * and Bottom() is the Vol() of a subbox split in the given direction
  195.  * and with the specified new upper bound.
  196.  */
  197.  
  198. long int Bottom(cube, dir, mmt)
  199. /* Compute part of Vol(cube, mmt) that doesn't depend on r1, g1, or b1 */
  200. /* (depending on dir) */
  201. struct box *cube;
  202. unsigned char dir;
  203. long int mmt[33][33][33];
  204. {
  205.     switch(dir){
  206.     case RED:
  207.         return( -mmt[cube->r0][cube->g1][cube->b1]
  208.             +mmt[cube->r0][cube->g1][cube->b0]
  209.             +mmt[cube->r0][cube->g0][cube->b1]
  210.             -mmt[cube->r0][cube->g0][cube->b0] );
  211.         break;
  212.     case GREEN:
  213.         return( -mmt[cube->r1][cube->g0][cube->b1]
  214.             +mmt[cube->r1][cube->g0][cube->b0]
  215.             +mmt[cube->r0][cube->g0][cube->b1]
  216.             -mmt[cube->r0][cube->g0][cube->b0] );
  217.         break;
  218.     case BLUE:
  219.         return( -mmt[cube->r1][cube->g1][cube->b0]
  220.             +mmt[cube->r1][cube->g0][cube->b0]
  221.             +mmt[cube->r0][cube->g1][cube->b0]
  222.             -mmt[cube->r0][cube->g0][cube->b0] );
  223.         break;
  224.     }
  225. }
  226.  
  227.  
  228. long int Top(cube, dir, pos, mmt)
  229. /* Compute remainder of Vol(cube, mmt), substituting pos for */
  230. /* r1, g1, or b1 (depending on dir) */
  231. struct box *cube;
  232. unsigned char dir;
  233. int   pos;
  234. long int mmt[33][33][33];
  235. {
  236.     switch(dir){
  237.     case RED:
  238.         return( mmt[pos][cube->g1][cube->b1] 
  239.            -mmt[pos][cube->g1][cube->b0]
  240.            -mmt[pos][cube->g0][cube->b1]
  241.            +mmt[pos][cube->g0][cube->b0] );
  242.         break;
  243.     case GREEN:
  244.         return( mmt[cube->r1][pos][cube->b1] 
  245.            -mmt[cube->r1][pos][cube->b0]
  246.            -mmt[cube->r0][pos][cube->b1]
  247.            +mmt[cube->r0][pos][cube->b0] );
  248.         break;
  249.     case BLUE:
  250.         return( mmt[cube->r1][cube->g1][pos]
  251.            -mmt[cube->r1][cube->g0][pos]
  252.            -mmt[cube->r0][cube->g1][pos]
  253.            +mmt[cube->r0][cube->g0][pos] );
  254.         break;
  255.     }
  256. }
  257.  
  258.  
  259. float Var(cube)
  260. /* Compute the weighted variance of a box */
  261. /* NB: as with the raw statistics, this is really the variance * size */
  262. struct box *cube;
  263. {
  264. float dr, dg, db, xx;
  265.  
  266.     dr = Vol(cube, mr); 
  267.     dg = Vol(cube, mg); 
  268.     db = Vol(cube, mb);
  269.     xx =  m2[cube->r1][cube->g1][cube->b1] 
  270.      -m2[cube->r1][cube->g1][cube->b0]
  271.      -m2[cube->r1][cube->g0][cube->b1]
  272.      +m2[cube->r1][cube->g0][cube->b0]
  273.      -m2[cube->r0][cube->g1][cube->b1]
  274.      +m2[cube->r0][cube->g1][cube->b0]
  275.      +m2[cube->r0][cube->g0][cube->b1]
  276.      -m2[cube->r0][cube->g0][cube->b0];
  277.  
  278.     return( xx - (dr*dr+dg*dg+db*db)/(float)Vol(cube,wt) );    
  279. }
  280.  
  281. /* We want to minimize the sum of the variances of two subboxes.
  282.  * The sum(c^2) terms can be ignored since their sum over both subboxes
  283.  * is the same (the sum for the whole box) no matter where we split.
  284.  * The remaining terms have a minus sign in the variance formula,
  285.  * so we drop the minus sign and MAXIMIZE the sum of the two terms.
  286.  */
  287.  
  288.  
  289. float Maximize(cube, dir, first, last, cut,
  290.         whole_r, whole_g, whole_b, whole_w)
  291. struct box *cube;
  292. unsigned char dir;
  293. int first, last, *cut;
  294. long int whole_r, whole_g, whole_b, whole_w;
  295. {
  296. register long int half_r, half_g, half_b, half_w;
  297. long int base_r, base_g, base_b, base_w;
  298. register int i;
  299. register float temp, max;
  300.  
  301.     base_r = Bottom(cube, dir, mr);
  302.     base_g = Bottom(cube, dir, mg);
  303.     base_b = Bottom(cube, dir, mb);
  304.     base_w = Bottom(cube, dir, wt);
  305.     max = 0.0;
  306.     *cut = -1;
  307.     for(i=first; i<last; ++i){
  308.     half_r = base_r + Top(cube, dir, i, mr);
  309.     half_g = base_g + Top(cube, dir, i, mg);
  310.     half_b = base_b + Top(cube, dir, i, mb);
  311.     half_w = base_w + Top(cube, dir, i, wt);
  312.         /* now half_x is sum over lower half of box, if split at i */
  313.         if (half_w == 0) {      /* subbox could be empty of pixels! */
  314.           continue;             /* never split into an empty box */
  315.     } else
  316.         temp = ((float)half_r*half_r + (float)half_g*half_g +
  317.                 (float)half_b*half_b)/half_w;
  318.  
  319.     half_r = whole_r - half_r;
  320.     half_g = whole_g - half_g;
  321.     half_b = whole_b - half_b;
  322.     half_w = whole_w - half_w;
  323.         if (half_w == 0) {      /* subbox could be empty of pixels! */
  324.           continue;             /* never split into an empty box */
  325.     } else
  326.         temp += ((float)half_r*half_r + (float)half_g*half_g +
  327.                  (float)half_b*half_b)/half_w;
  328.  
  329.         if (temp > max) {max=temp; *cut=i;}
  330.     }
  331.     return(max);
  332. }
  333.  
  334. int
  335. Cut(set1, set2)
  336. struct box *set1, *set2;
  337. {
  338. unsigned char dir;
  339. int cutr, cutg, cutb;
  340. float maxr, maxg, maxb;
  341. long int whole_r, whole_g, whole_b, whole_w;
  342.  
  343.     whole_r = Vol(set1, mr);
  344.     whole_g = Vol(set1, mg);
  345.     whole_b = Vol(set1, mb);
  346.     whole_w = Vol(set1, wt);
  347.  
  348.     maxr = Maximize(set1, RED, set1->r0+1, set1->r1, &cutr,
  349.             whole_r, whole_g, whole_b, whole_w);
  350.     maxg = Maximize(set1, GREEN, set1->g0+1, set1->g1, &cutg,
  351.             whole_r, whole_g, whole_b, whole_w);
  352.     maxb = Maximize(set1, BLUE, set1->b0+1, set1->b1, &cutb,
  353.             whole_r, whole_g, whole_b, whole_w);
  354.  
  355.     if( (maxr>=maxg)&&(maxr>=maxb) ) {
  356.     dir = RED;
  357.     if (cutr < 0) return 0; /* can't split the box */
  358.     }
  359.     else
  360.     if( (maxg>=maxr)&&(maxg>=maxb) ) 
  361.     dir = GREEN;
  362.     else
  363.     dir = BLUE; 
  364.  
  365.     set2->r1 = set1->r1;
  366.     set2->g1 = set1->g1;
  367.     set2->b1 = set1->b1;
  368.  
  369.     switch (dir){
  370.     case RED:
  371.         set2->r0 = set1->r1 = cutr;
  372.         set2->g0 = set1->g0;
  373.         set2->b0 = set1->b0;
  374.         break;
  375.     case GREEN:
  376.         set2->g0 = set1->g1 = cutg;
  377.         set2->r0 = set1->r0;
  378.         set2->b0 = set1->b0;
  379.         break;
  380.     case BLUE:
  381.         set2->b0 = set1->b1 = cutb;
  382.         set2->r0 = set1->r0;
  383.         set2->g0 = set1->g0;
  384.         break;
  385.     }
  386.     set1->vol=(set1->r1-set1->r0)*(set1->g1-set1->g0)*(set1->b1-set1->b0);
  387.     set2->vol=(set2->r1-set2->r0)*(set2->g1-set2->g0)*(set2->b1-set2->b0);
  388.     return 1;
  389. }
  390.  
  391.  
  392. Mark(cube, label, tag)
  393. struct box *cube;
  394. int label;
  395. unsigned char *tag;
  396. {
  397. register int r, g, b;
  398.  
  399.     for(r=cube->r0+1; r<=cube->r1; ++r)
  400.        for(g=cube->g0+1; g<=cube->g1; ++g)
  401.       for(b=cube->b0+1; b<=cube->b1; ++b)
  402.         tag[(r<<10) + (r<<6) + r + (g<<5) + g + b] = label;
  403. }
  404.  
  405. doquant()
  406. {
  407. struct box    cube[MAXCOLOR];
  408. unsigned char    *tag;
  409. int        next;
  410. register long int    i, weight;
  411. register int    k;
  412. float        vv[MAXCOLOR], temp;
  413.  
  414.     /* input R,G,B components into Ir, Ig, Ib;
  415.        set quant_size to width*height */
  416.         /* set K to number of colors in quantized image */
  417.  
  418.   bzero (m2, 33*33*33*sizeof (float));
  419.   bzero (wt, 33*33*33*sizeof (long int));
  420.   bzero (mr, 33*33*33*sizeof (long int));
  421.   bzero (mg, 33*33*33*sizeof (long int));
  422.   bzero (mb, 33*33*33*sizeof (long int));
  423.  
  424.     Hist3d(wt, mr, mg, mb, m2);
  425.     free(Ig); free(Ib); free(Ir);
  426.  
  427.     M3d(wt, mr, mg, mb, m2);
  428.  
  429.     cube[0].r0 = cube[0].g0 = cube[0].b0 = 0;
  430.     cube[0].r1 = cube[0].g1 = cube[0].b1 = 32;
  431.     next = 0;
  432.         for(i=1; i<K; ++i){
  433.             if (Cut(&cube[next], &cube[i])) {
  434.               /* volume test ensures we won't try to cut one-cell box */
  435.               vv[next] = (cube[next].vol>1) ? Var(&cube[next]) : 0.0;
  436.               vv[i] = (cube[i].vol>1) ? Var(&cube[i]) : 0.0;
  437.         } else {
  438.               vv[next] = 0.0;   /* don't try to split this box again */
  439.               i--;              /* didn't create box i */
  440.         }
  441.             next = 0; temp = vv[0];
  442.             for(k=1; k<=i; ++k)
  443.                 if (vv[k] > temp) {
  444.                     temp = vv[k]; next = k;
  445.         }
  446.             if (temp <= 0.0) {
  447.               K = i+1;
  448.               /* fprintf(stderr, "Only got %d boxes\n", K); */
  449.               break;
  450.         }
  451.     }
  452.  
  453.     /* the space for array m2 can be freed now */
  454.  
  455.     tag = (unsigned char *)malloc(33*33*33);
  456.     for(k=0; k<K; ++k){
  457.         Mark(&cube[k], k, tag);
  458.         weight = Vol(&cube[k], wt);
  459.         if (weight) {
  460.         lut_r[k] = Vol(&cube[k], mr) / weight;
  461.         lut_g[k] = Vol(&cube[k], mg) / weight;
  462.         lut_b[k] = Vol(&cube[k], mb) / weight;
  463.         }
  464.         else{
  465.           /* fprintf(stderr, "bogus box %d\n", k); */
  466.           lut_r[k] = lut_g[k] = lut_b[k] = 0;        
  467.         }
  468.     }
  469.  
  470.     for(i=0; i<quant_size; ++i) Qadd[i] = tag[Qadd[i]];
  471.  
  472.         free (tag);
  473.  
  474.         return 0;
  475.     
  476.     /* output lut_r, lut_g, lut_b as color look-up table contents,
  477.        Qadd as the quantized image (array of table addresses). */
  478. }
  479.